Localized immune surveillance of primary melanoma in the skin deciphered through executable modeling

While skin is a site of active immune surveillance, primary melanomas often escape detection. Here, we have developed an in silico model to determine the local cross-talk between melanomas and Langerhans cells (LCs), the primary antigen-presenting cells at the site of melanoma development. The model predicts that melanomas fail to activate LC migration to lymph nodes until tumors reach a critical size, which is determined by a positive TNF-α feedback loop within melanomas, in line with our observations of murine tumors. In silico drug screening, supported by subsequent experimental testing, shows that treatment of primary tumors with MAPK pathway inhibitors may further prevent LC migration. In addition, our in silico model predicts treatment combinations that bypass LC dysfunction. In conclusion, our combined approach of in silico and in vivo studies suggests a molecular mechanism that explains how early melanomas develop under the radar of immune surveillance by LC.


Langerhans cell model description
The Langerhans Cell (LC) Model is designed to describe how LCs respond to various cytokines and growth factors. The outputs of the model are Residency, Proliferation and Survival. All nodes in the LC network vary between 0-2. As LCs migrate, proliferate and die at a given rate in the epidermis (97), each of the output nodes will have value 1 in the unperturbed state. An increase to level 2 in Residency corresponds to a loss of normal migration while a decrease to level 0 indicates enhanced migration. Note that, where two nodes with the same name (such as MEK) exist in both the LC and melanoma models, the suffix "_LC" has been added to the LC version of the node.
TGF-β is a key cytokine controlling LC behaviour. Loss of TGF-β in the epidermis causes migration of LCs out of the skin, leaving the epidermis devoid of LCs (62). TGF-β signalling is mediated by the TGF-β receptor which mainly signals via the Smad2/3 complex and LAMTOR-p14 in LCs. Smad signalling engages activation of Pu.1, RUNX3 and Id2 transcription factors, which is key to maintain residency in the epidermis (98). TGF-β also activates the LAMTOR-p14 complex to maintain ERK and mTOR activity in LCs (99) for proliferation and survival signalling. Finally, TGF-β signalling contributes to expression of E-cadherin and other adhesion molecules such as EpCAM and Claudin-1 (100,101).
LCs express the CSF1R, Axl and Tyro3 receptor tyrosine kinases, which all contribute to proliferative and survival signalling through PI3K. CSF1R is a receptor both for CSF1 and IL-34; IL-34 provides vital survival signals for LCs in the epidermis (102). Axl is expressed downstream of TGF-β signalling and feedback between Axl and Tyro3 means that Tyro3 is expressed only in the absence of Axl (103). The ligands for Tyro3 are Gas6 and Pros1, while Axl accepts only Gas6 (100).
LCs respond to Wnt ligands through the canonical Wnt pathway involving the Frizzled receptor and activation of β-catenin (104,105), contributing to proliferative signalling. β-catenin is also regulated through interaction with E-cadherin (106), meaning that E-cadherin down-regulation by TGF-β signalling results in β-catenin activation (101). DKK proteins are soluble inhibitors of WNT-Frizzled interactions (105) and LCs are also responsive to DKK (104).
TNF-α is a key danger signal in the epidermis and LCs respond mainly through the TNFRII receptor (57). Keratinocytes produce the inflammatory cytokines IL-1β and TNF-α (in response to IL-1β) (57). These two signals induce LC migration (107) and converge on the IRF1 transcription factor (108). Co-expression of IRF1 with IRF4 is also associated with proinflammatory migratory LCs (11). LC migration requires downregulation of E-cadherin (109).

Melanoma model description
The core melanoma model aims to be a representation of melanoma cells that accounts for the main drivers of cutaneous melanoma and can use these drivers to predict their response (in terms of proliferation and apoptosis) to a variety of external stimuli and targeted therapies. It focuses on BRAF-and NRAS-driven tumours as these make up a significant fraction of cutaneous melanomas (43) and are the focus of nearly all experimental investigations. While the model has not been tested on NF1-driven tumours, due to a lack of representation in literature, it is likely it would behave well as NF1 is an inhibitor of RAS (110), and NRAS-driven tumours have been tested. "Triple-wild-type" tumours without mutations in any of these genes make up a small but significant fraction of melanomas but are thought to be a heterogeneous group (43), and it is unclear how well the model would be able to predict the behaviour of any of these melanomas.
Melanoma cell lines or pre-clinical models are simulated using the known driver mutations present in these cells (Supplementary Table 13). Note that while melanomas do not always show CDKN2A/B mutations, they nearly always show either mutations or loss of expression of these genes (77). Therefore, we set the CDKN2AB node to 0 for all melanoma cells. The model attempts to simulate the impact of MITF (varying between 0 and 3) on the transcriptional state (44) but regulation of the MITF state itself is considered outside of the scope of the model. Information on the MITF-status of cell lines or tumours is not always provided so it is assumed that melanoma cells are in the MITF high (MITF = 2) state unless there is information to the contrary.
Each of the nodes in the network is assigned to one of four "pathways" generally describing their behaviour, we will briefly outline each of these below.

Growth Signalling pathway
The growth pathway describes general growth signalling, usually mediated by growth factors and receptor tyrosine kinases, but commonly aberrantly activated in melanoma. This consists of the MAPK and PI3K pathways and their downstream targets. The model includes two inputs which are necessary for growth signalling -a generic growth factor and cell-cell contact, which is required for the maintenance of normal melanocytes. This growth factor is kept generic as many factors have been implicated either in the normal growth of melanocytes (for example ET-1 and FGF (111)) or in resistance to targeted therapies (for example by upregulation of RTKs like EGFR or PDGFR (112)), and they all converge on the same pathways (113). The level of growth factors ranges from 0-4, accounting for starvation (0), levels in standard media (1), levels used to stimulate melanocyte growth (2) and high levels that can drive melanoma growth (4). Activation of RTKs and FAK (via cell-cell contact signalling) leads to activation of PI3K and RAS, which activates both PI3K and RAF.
The MAPK pathway progresses by activation of RAF, MEK and ERK. While feedback from the MAPK pathway has been described (114), we have been unable to link these effects to the changes in steady state behaviour in terms of apoptosis and proliferation that are functionally important for this kind of modelling and so this has been excluded. In order to represent the complex interaction between different mutations observed in melanomas with acquired resistance to targeted therapies, MAPK pathway mutations are handled differently to other mutations in the model. Generally, a gain-of-function mutation will mean we set the value of a specific node to its highest value and to model drugs we set the value to 0. For MAPK node we include additional nodes to represent mutations and drugs, which then influence the node itself. For example, in addition to a "BRAF" node, we have a "BRAFV600E" node representing the gain-of-function mutation, a "BRAFamplification" node representing copy number gain of the mutant BRAF and a "BRAFi" node representing a mutant BRAF inhibitor like dabrafenib. The "BRAFV600E" node induces RAS-independent BRAF activity, while the "BRAFamplification" node enhances the degree of RAS-independent BRAF activity. The "BRAFi" node specifically inhibits RAS-independent BRAF activation. This applies to MEK and ERK but not NRAS.
The PI3K pathway includes PTEN, a tumour suppressor gene frequently mutated in melanoma, and an "Akt" node which represent any of the Akt isoforms. Downstream of the PI3K pathway lies mTORC1 and 2, S6K and GSK3b. The Growth Signalling pathway also includes Myc which is downstream of many pathways including ERK, β-catenin, GSK3b, TGF-β and p53.
The Growth Signalling pathway is largely stand-alone, but it does have some feedback from the rest of the model. The RTK node can be upregulated by cJUN, to model the known upregulation of RTKs (42,112,115) in the MITF-low transcriptional state, which is characterised by high cJUN activity (48). PTEN is regulated by p53 and cJUN (116) and RAF is activated by PKC to model the known action of the tumour promoter TPA (12-O-Tetradecanoylphorbol-13-acetate (117)).

Intercellular Signalling pathway
The secreted factors pathway describes how signals from TGF-β, TNF-α, α-MSH and WNT ligands are recognised by melanoma cells and how these signals are propagated. Some of these factors, such as WNT ligands and TGF-β are produced by melanoma cells. TGF-β expression is downstream of the MAPK pathway (118) and so is controlled by ERK in the model. The control of expression of WNT ligands is less well understood and so, except for WNT5a which is inhibited by MITF to mirror its association with the low MITF state (119), the WNT ligands are inputs to the model. TNF-α can be produced in small amounts by melanoma cells in culture (36), but has been thought to be primarily derived from other elements of the tumour environment. However, in the context of the epidermis, we found that melanoma cells had the highest level of TNF expression ( Figure 3D). The original melanoma model allowed for either possibility by treating the level of TNF-α as an input to the model, while the melanoma-LC includes an in-depth treatment of TNF expression (detail provided in main text). α-MSH is produced by keratinocytes (120) and is sometimes included in media to culture melanocytes and so is considered an input to the model.
In the base melanoma model, TNF-α signals via the TNFR node and stimulates IKK, leading to activation of NF-κB, although activation of IKK through Akt and ERK and synergistic activation of NF-κB targets by cJUN are also included in the model (48). Although TNF-α is known to induce the MITF-low state (48,121), the exact mechanism behind this switch is poorly understood and in any case, mechanisms of phenotype switching are considered beyond the scope of the model, and so this effect must be taken into account at the point of specifying inputs. TNF-α is known to induce apoptosis in some cell types however it appears to have the opposite effect in melanoma, and in some contexts it even protects from apoptosis (122). TNF signalling also activates JNK and p38 (123).
The TGF-β pathway is modelled by its receptor, TBR (TGF-βR), and Smad2/3, whose effects are broadly anti-proliferative in melanoma. The receptor TBR is inhibited by MITF to account for the fact that MITF low cells have enhanced TGF-β signalling (124) and induction of low MITF state causes induction of TGFBR2 gene (112), while conversely MITF overexpression inhibits TGF-β signal propagation (125).
The α-MSH pathway signals through the MC1R receptor leading to production of cAMP, and activation of PKA and CREB (126). CREB can also be activated by high levels of ERK activity (127) and its downstream target in the model is cJUN (128). PKA is also negatively regulated by Smad2/3 (129) and can activate β-catenin (130).
In melanoma, both canonical and non-canonical WNT signalling have been reported (119). Canonical WNT signalling through the Frizzled receptor leads to β-catenin activation through inhibition of Axin, while multiple canonical WNT ligands (WNT3a, WNT10a, etc) exist they are not generally distinguished in the literature and so a single WNT node is used for canonical WNT ligands. The non-canonical WNT ligand Wnt5a signals through the ROR2 receptor, stimulating PKC and inhibiting the canonical WNT pathway through activation of the Siah2 E3 ligase which targets β-catenin for destruction (131). Generally, in melanoma, the two forms of WNT signalling are found to be mutually exclusive and in particular associated with transcriptional state (124). Wnt5a in particular is associated with and can induce the MITF low state (132), and therefore in the model is inhibited by MITF. Little appears to be known about regulation of WNT ligand expression in melanoma, but we have included the DKK proteins which are intercellular inhibitors of WNT receptor binding that are known to be differentially expressed in melanomas (39). The role of β-catenin is not totally clear in melanoma, with some studies arguing for a pro-melanoma role (133,134) while others argue that loss of β-catenin can enhance melanoma growth (135). These contradictions are likely due to the interplay between β-catenin activity and transcriptional state (136), where MITF itself has also been seen to have both pro-and anti-proliferative effects in melanoma (137,138). However, as regulation of MITF is considered beyond the scope of the model, we have not attempted to model the effects of β-catenin in melanoma progression and the only downstream targets of β-catenin, through TCF/LEF1 are p16 and Myc.

Cell Cycle pathway
The cell cycle pathway controls proliferation by modelling the G1/S checkpoint as a function of Cyclin D and CDK4 activity. While other cyclins, CDKs and checkpoints exist, the evidence suggests that this checkpoint is usually the limiting factor in proliferation of melanoma cells. Proliferation is downstream of E2F, which is controlled primarily by pRB which in turn is regulated by the CycD-CDK4 complex. Note that both the proliferation and apoptosis nodes very between 0 and 4 should be interpreted at the population level, so a higher level of proliferation may suggest either a higher fraction of proliferating cells or a shorter time between cell division and similarly for apoptosis.
The tumour suppressor genes p14/ARF (p19ARF in mice), p15 and p16/INK4A are included in the model with p15 and p16 inhibiting CDK4 and ARF targeting both Mdm, which regulates p53, and E2F (139). These proteins are expressed from the CDKN2A and CDKN2B genes which are close together in the genome and are frequently mutated in melanoma (43,118). Many early papers debated whether genetic loss of factor such as p16 was critical for melanoma development (140,141). However, it has become increasingly clear that culturing of melanocytes can select for loss of p16 expression at the transcriptional level without underlying mutations (142). Furthermore, when a panel of uncultured melanoma cells were examined using DNA and RNA profiling, nearly all had either mutated CDKN2A, lost expression of p16 or had an additional known compensatory mutation in the CDK4 pathway (e.g. CDK4R24C, (77)). For these reasons, we assume that loss of CDKN2AB at either the genetic or transcriptional level is a pre-requisite for melanoma formation, and when simulating cell lines, we assume a loss of CDKN2AB even when none is recorded in the literature (as this may be at the transcriptional, not genetic level).
The tumour suppressors p21 and p27 are also included and their activity converges on the CycD-CDK4 complex. p21 lies downstream of p53, MITF, FRA1 and JUNB and is inhibited by Akt and TBX2. p53 itself is regulated by its inhibitor Mdm and a node called "ReplicationStress", which is used to distinguish normal from transformed melanocytes. Although the exact details of this process are not well understood, this distinction is necessary to understand for example the differential sensitivity of melanocytes and melanoma cells to BH3 mimetic drugs (143). Apart from p21, p53 also regulates PUMA and BaxBak, in the survival pathway. p27 is assumed to have a level of constitutive activity in the model but is also activated by Smad2/3 and is inhibited by MITF, Akt, Myc and p90RSK.
MITF is also included in the cell cycle pathway and is an input to the model, as its regulation is highly complex (144) and for now we consider it outside the scope of the model. MITF is considered to act as a rheostat in melanoma, having quite different effects at its different levels of activity. At a very high level of MITF expression, cells exhibit a highly differentiated phenotype that is resistant to cell death but is also slowly proliferative. At a high level of MITF expression (which can be considered the "normal" level for melanoma cells), cells are proliferative. At a low level of MITF expression, cells become de-differentiated and are more mobile and less proliferative but more resistant to therapy due to upregulation of RTKs. Complete loss of MITF expression is highly detrimental to cells and is not seen in melanoma cells without experimental perturbation. It is worth noting that cells can switch between the MITF-high and MITF-low without mutation and that any population of cells will likely include some of each type (145). Furthermore, recent single-cell RNAseq experiments have dissected these states in even greater detail, and in fact it is possible that further subtypes dependent on other factors such as SOX10 may define the phenotype in greater detail (55,(146)(147)(148)(149). However, in the interests of keeping the model focused, we decided to limit the model's treatment of MITF to a simple rheostat. MITF regulates many other proteins in the model, which will be described as the arise in the text. One key downstream target of MITF is cJUN, which is thought to positively define the transcriptional program underlying the MITF low state (48,150). cJUN has various pro-melanoma roles, including upregulation of RTKs, inhibition of PTEN and co-activation of NF-κB targets. Note that the model does not deal with senescence in melanoma. Short-term growth followed by stable growth arrest in naevi, which usually remain benign for decades, is a characteristic feature of melanocytes with oncogenic mutations (151). However, studies into senescence are plagued by irreproducibility and inconsistent definitions, at least in part due to significant differences between in vivo and in vitro conditions. As tools have developed to a point where naevi can be studied in vivo, a picture has emerged suggesting naevi formation is mediated by secreted factors rather than cell intrinsic factors (24). While p14, p15 and p16 are upregulated in naevi, there is significant evidence that genetic or transcriptional changes beyond loss of these factors are required for transition from naevus to melanoma (141,152,153). Longitudinal profiling has helped to identify CDKN2A and TERT promoter mutations as unlikely to occur in naevi but common in melanoma (77). For now, we largely follow McNeal et al. (2015) (118) for our treatment of naevi, but this is not a focus for the model and so they are largely avoided in the specification. Note that we also avoid discussion of TERT promoter mutations as these are complex to model, they are only important after tens of rounds of divisions, and they simply allow for long term growth without interacting with other parts of the melanoma network.

Survival pathway
The survival pathway describes how signals from the rest of the model are translated into apoptotic signals in the cell. As described above, the primary pro-apoptotic signal in the cell is activation of BaxBak by p53, which leads to mitochondrial outer membrane permeabilization (MOMP). This process is opposed by Bcl-2, Mcl-1 and BCL2A1, which are inhibited by PUMA, BIM and BAD, which are, in turn inhibited by growth signalling from the MAPK and PI3K pathways.
MOMP leads to cell death by both caspase-dependent and caspase-independent pathways. The caspase-dependent pathway is initiated by caspase 3 activation but in melanoma is often not significant (154) due to expression of the caspase inhibitor ML-IAP (155). In order to balance loss of ML-IAP expression in the MITF low state, we have included an edge linking NF-κB to expression of IAP (a node comprising cIAP1/2), a known NF-κB target in fibroblasts (156). Caspase-independent death is mediated by AIF, which can also be activated by extremely high levels of MAPK signalling via PARP and PAR (157).

Melanoma model specification deviations
Line numbers correspond to Supplementary Table S6. Note that "numerical deviations", where the overall behaviour of the system is maintained but there are differences in the exact level of the node, are not described in detail here.  (118) state that p16 overexpression has no effect on growth of melanoma cells, in the model this has an effect. This is because the role of p16 in melanoma is unclear, but expression is frequently lost, suggesting a functional role.

Melanoma-LC model specification deviations
Line numbers correspond to Supplementary Table S9. These issues are in addition to those outlined above, unless otherwise indicated. Unlike those listed above, these issues are "numerical deviations" where the behaviour of the system is not seriously affected, but the exact numbers do not match.
• Lines 43, 45, 47, 49: Zhuang et al, (2008) (158) found that knockdown of Myc with shRNA leads to reduction of proliferation in both the SK-MEL-2 and SK-MEL-19 backgrounds, with the same effect occurring when combined with p53 shRNA. The melanoma-LC model finds two steady states when simulating these experiments, in one the correct behaviour is identified and in the other the TNF-α signalling loop keeps JNK activity high, leading to enhanced cJUN activity and reducing the impact of Myc shRNA. As these experiments were performed in cell culture, and we found that YUMM1.7 cells in cell culture do not produce substantial TNF-α, it is possible that their measurements were of cells without engagement of the signalling loop. • Line 71: Posch et al, (2013) (161) found that a combined PI3K/mTORC1/2 inhibitor (GSK2126458) reduced proliferation of d04 cells. Similar to above, our simulations of the melanoma-LC model found two steady states divided on engagement of the TNF-α signalling loop. In the state with TNF-α signalling engaged, a low level of proliferation is maintained.

Supplementary Data
SupplementaryData1.xlsx Fixed points of the network under single and double perturbations of the network, as determined by computational screening. Perturbations refer either to alterations to specific node (level specified in data file) or simulation of targeted therapy (levels of affected nodes specified in Supplementary Table S11). Each line refers to the effect of a specific perturbation on a single output node (melanoma cell apoptosis or proliferation or Langerhans cell apoptosis, proliferation or residency).      Supplementary Table S11.  Caspase-dependent apoptosis can be initiated via intrinsic pathway and is inhibited by IAPs.           Supplementary Table S7. Regulatory edges in the melanoma-LC model. Each row shows a single edge, originating from the node in the "from" column and directed to the node in the "to" column. The table also contains details of the type of edge, a brief explanation of the mechanism of this regulation and a reference to the literature where required.    Frizzled is receptor for Wnt, can be inhibited by DKK, but at endogenous levels Wnt signalling still occurs.